class BlockEP
{
	public:
		BlockEP(mat L, mat b)
		{
			_L=L;
			_p=L.n_cols;
			d=new boost::math::normal_distribution<>(0,1); 
			boost::random::uniform_01<> U;
			_runif=new RandomG::Random<boost::random::uniform_01<> >(U);
			_b=b;
			mat t(2,1);
			G=t;
			mat I(2,2);
			I.fill(1);
			_I=diagmat(I);
		}
		void operator()(int j, mat x,int k, double *z, double *w)
		{
			double tr1=_b(0,j);
			double tr2=_b(0,j+1);
			for(int i=0;i<j;i++) 
			{
				tr1-=_L(j,i)*x(0,i);
				tr2-=_L(j+1,i)*x(0,i);
			}
			tr1=tr1/_L(j,j);
			if(k==0)
			{
				G(0,0)=_L(j+1,j);
				G(1,0)=_L(j+1,j+1);
				norm=sqrt(dot(G,G));
				norm3=norm*norm*norm;
			}
			double cd= cdf(*d,tr2/norm);
			if(cd > .999999999) cd=.999999999;
			if(cd < .0000000001) cd=.0000000001;
			mat dlZm=((double)1/(1-cd)*pdf(*d,tr2/norm))*G/norm;
			mat dlZS=((double)1/(1-cd)*pdf(*d,tr2/norm)*tr2/(norm3))*G*G.t();


			mat Sig=_I-dlZm*dlZm.t()+2*dlZS;
			double sig=sqrt(Sig(0,0));
			double m=dlZm(0,0);	
			double cd2=cdf(*d,(tr1-m)/sig);	
			double u = (*_runif)(); 
			double arg=u+(1.-u)*cd2;
			if(arg > .999999999) arg=.999999999;
			if(arg < .0000000001) arg=.0000000001;
			*z = quantile(*d,arg)*sig+m;
			*w=(1-cd2)*pdf(*d,*z)*sig/(pdf(*d,(*z-m)/sig));

		}
	private:
		mat _L;
		int _p;
		mat _b;
		mat _I;
		boost::math::normal_distribution<> *d;
		RandomG::Random<boost::random::uniform_01<> > *_runif;
		mat G;
		double norm3;
		double norm;
};


